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Molecular dynamics (MD) simulations of ions (K+, Na^, Ca^+ and Cl“) in aqueous 
solutions are investigated. Water is described using the SPC/E model. A stochastic 
coarse-grained description for ion behaviour is presented and parameterized using MD 
simulations. It is given as a system of coupled stochastic and ordinary differential 
equations, describing the ion position, velocity and acceleration. The stochastic coarse¬ 
grained model provides an intermediate description between all-atom MD simulations 
and Brownian dynamics (BD) models. It is used to develop a multiscale method which 
uses all-atom MD simulations in parts of the computational domain and (less detailed) 
BD simulations in the remainder of the domain. 
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1. Introduction 

Molecular dynamics (MD) simulations of ions in aqueous solutions are limited 
to modelling processes in relatively small domains containing (only) several 
thousands of water molecules ai. Ions play important physiological functions 
in living cells which typically consist of lO^^-IO^^ water molecules. In particular, 
processes which include transport of ions between different parts of a cell cannot 
be simulated using standard all-atom MD approaches. Coarser models are instead 
used in applications. Examples include Brownian ^namics (BD) simulations Q 
and mean-field Poisson-Nernst-Planck equations y|. In BD methods, individual 
trajectories of ions are described using 

dXi = VwdWi, i = l,2,3, (1.1) 

where X = \Xi,X2,X^ is the position of the ion, D is its diffusion constant and 
Wi, i = 1,2,3, are three independent Wiener processes jb]. BD description (II.ip 
does not explicitly include solvent molecules in the simulation. Moreover, 
in applications, equation dni can be discretized using a (nanosecond) time 
step which is much larger than the typical time step of MD simulations 
(femtosecond) Q. This makes BD less computationally intensive than the 
corresponding MD simulations. 
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Longer time steps of BD simulations enable efficient simulations of ion 
transport between different parts of the cell, but they limit the level of detail 
which can be incorporated into the model. For example, intracellular calcium 
is regulated by the release of Ca^'*' ions from the endoplasmic reticulum via 
inisitol-4,5-triphosphate receptor (IP 3 R) channels. BD models in the literature use 
equation (II.ip to describe trajectories of calcium ions i 0 . The conformational 
changes between the open and closed states of IP 3 R channels are controlled by the 
binding of Ca^'*' to activating and inhibitory binding sites. BD models postulate 
that binding of an ion occurs with some probability whenever the distance between 
the ion and an empty site is less than the specific distance, the so called reaction 
radius i, i- Although details of the binding process are known 113, in, they 
cannot be incorporated into coarse BD models of calcium dynamics, because 
equation dni) does not correctly describe short time behaviour of ion dynamics. 

The calcium induced calcium release through IP 3 R channels is an example 
of a multiscale dynamical problem where MD simulations are important only in 
certain parts of the computational domain (close to an IP 3 R channel), whilst in 
the remainder of the domain a coarser, less detailed, BD method could be used (to 
describe trajectories of ions). Such multiscale problems cannot be simulated using 
MD methods, but there is potential to design multiscale computational methods 
which compute the desired information with an MD-level of resolution by using 
MD and BD models in different parts of the computational domain (l3 |. 

In jl2|, three relatively simple and analytically tractable MD models are 
studied (describing heat bath molecules as point particles) with the aim of 
developing and analyzing multiscale methods which use MD simulations in parts 
of the computational domain and less detailed BD simulations in the remainder 
of the domain. In this follow up paper, the same question is investigated in all¬ 
atom MD simulations which use the SPC/E model of water molecules. In order to 
couple MD and BD simulations, we need to first show that the MD model is in a 
suitable limit described by a stochastic model which does not explicitly take into 
account heat bath (water) molecules. In jl^ . this coarser description was given in 
terms of Langevin dynamics. Considering all-atom MD simulations, the coarser 
stochastic model of an ion is more complicated than Langevin dynamics. In this 
paper, it will be given by 


dAi 

= Vidt, 


( 1 . 2 ) 

dCi 

= Uidt, 


(1.3) 

dUi 

= {-myi + Zi)dt, 


(1.4) 

dZi 

= -imZi + ri3Ui)dt + r]4dWi, 

i — 1,2, 3, 

(1.5) 


where X= [Ai, A 2 ,X 3 ] is the position of the ion, V = [Vi,I^, V 3 ] is its velocity, 
U = [C/i, 1 / 2 , is its acceleration, Z = [Zi, Z 2 , .^ 3 ] is an auxiliary variable, 
dW = [dITi, dlL 2 , dIT 3 ] is white noise and rjj, j = 1, 2,3,4, are parameters. These 
parameters will be chosen according to all-atom MD simulations as discussed in 
Section |3l In Section HJ we show that (|1.2I) - (|1.5I) provides a good approximation 
of ion behaviour. In Section [5l we further analyse the system (ll.2p - (ll.5l) and show 
how parameters r/j, j = 1, 2, 3,4, can be connected with diffusion constant D used 
in the BD model (|l.ll) . 
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ion 

Aio 

[Da A^^ ps~^] 

Bio 

[Da A® ps“^] 

qi 

N 

M 

[Da] 

K+ 

2.8973 X 10® 

2.4587 X 10^ 

+1 

39.0983 

Na+ 

6.6813 X lO’^ 

1.1807 X 10^ 

+1 

22.9898 

Ca2+ 

1.1961 X 10® 

1.5797 X 10^ 

+2 

40.078 

cr 

1.8038 X 10^ 

6.1347 X 10^ 

-1 

35.453 


Table 1. Parameters of all-atom MD simulations of ions. 


The coarse-grained model (|1.2I) - (|1.5I) is used as an intermediate model between 
the all-atom MD model and BD description dni). In Section |5l we show how it can 
be coupled with the BD model which uses a much larger time step than the MD 
model. In Section |6l the coarse-grained model is coupled with all-atom 

MD simulations. We then show that all-atom MD models of ions can be coupled 
with BD description (HU) using the intermediate coarse-grained model m- 
(HU- We conclude by discussing related methods developed in the literature in 
Section [3 


2. Molecular dynamics simulations of ions in SPC/E water 


There have been several MD models of liquid water developed in the literature. 
The simplest models (for example, SPC [T^, SPC/E jl4| and TIP3P [IH) include 
three sites in total, two hydrogen atoms and an oxygen atom. More complicated 
water models include four, five or six sites EES- In this paper, we use the three- 
site SPC/E model of water which was previously used for MD simulations of ions 
in aqueous solutions (T^ . [l|. In the SPC/E model, the charges (g/i = 0.4238e) 
on hydrogen sites are at lA from the Lennard-Jones center at the oxygen site 
which has negative c harg e qo = —0.8476e. The HOH angle is 109.47°. We use the 
RATTLE algorithm [l(^ to satisfy constraints between atoms of the same water 
molecule. 

We investigate four ions (K'*', Na"*", Ca^"*" and Cl“) at 25°C using MD 
parameters presented in [l^ . Let us consider a water molecule and let us denote by 
Tjo (resp., rn and rj 2 ) the distance between the ion and the oxygen site (resp., the 
first and second hydrogen sites). The pair potential between the water molecule 
and the ion is then given by P, , 


A,- 




+ ke 


qm 

Til 


+ ke 


qm 

D2 ’ 


( 2 . 1 ) 


where Ajo and Bio are Lennard-Jones parameters between the oxygen on the 
water molecule and the ion, ke is Coulomb’s constant and Qi is the charge on the 
ion. The values of parameters are given for four ions considered in Table [T] We 
express mass in daltons (Da), length in angstroms (A) and time in picoseconds 
(ps), consistently in the whole paper. Using these units, the parameters of the 
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ion 

D 

[A^ ps“^] 

{vh 

[A^ ps“^] 

m 

[A^ ps“^] 

(Z?) 

|A2 ps-«| 

K+ 

0.183 

6.32 

4.86 X 10® 

1.65 X 10^ 

Na+ 

0.128 

10.8 

2.21 X 10^ 

8.88 X 10^ 

Ca2+ 

0.053 

6.18 

1.87 X 10"^ 

9.23 X 10^ 

cr 

0.177 

6.98 

6.56 X 10® 

2.97 X 10^ 


Table 2. Average values obtained by all-atom MD simulations. 


Lennard-Jones potential between the oxygen sites on two SPC/E water molecules 
are Aqq = 2.6334 x 10^ Da ps“^ and Bqo = 2.6171 x 10® Da A® ps“^. 

We consider a cube of side L = 24.83 A containing 511 water molecules and 
1 ion, i.e. we have 8® = 512 molecules in our simulation box. In the following 
section, we use standard NVT simulations where the temperature is controlled 
using Nose-Hoover thermostat (^.1^ and the number of particles is kept constant 
by implementing periodic boundary conditions. In particular, we assume that our 
simulation box is surrounded by periodic copies of itself. Then the long-range 
(Coulombic) interactions can be computed using several different approaches, 
including the Ewald summation or the reaction held method j^, [^. We use 
the cutoff sphere of radius L/2 and the reaction held correction as implemented 
in [2|. This approach is more suitable for multiscale methods (studied later in 
Section [6| than the Ewald summation technique. The MD timestep is for all MD 
simulations in this paper chosen as At = 10“® ps = 1 fs. 


3. Parametrization of the coarse-grained model of ion 

In MD simulations, an ion is descibed by its position X = [Xi, X2, X^] and velocity 
V = [Vi,V 2 , V 3 ] which evolve according to 

dXi = Vidt, (3.1) 

MdVi = Fidt, i = 1,2,3, (3.2) 

where M is the mass of the ion (given in Table [T]) and F = [Fi, F2, F3] is the 
force acting on the ion. We use all-atom MD simulations as described in Section [2] 
to estimate diffusion coefficient D and second moments of V) and Ui = Fi/M, 
i = l,2, 3. They are given for four ions considered in Table[2] To estimate {Uf), we 
calculate the average force in the Tth direction (F^) where (•) denotes an average 
over sufficiently large time interval (nanosecond) of MD simulations. Taking into 
account the symmetry of the problem, we estimate (Uf) = {Ff) /M^ as the average 
over all three dimensions 

{UD + m) + {UD 

3 

This value is reported in Table |2j In the same way, the reported values of {V^) 
are computed as averages over all three dimensions. Diffusion constant D can be 
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ion 

Vi 

[ps-2] 

m 

[ps-i] 

V3 

[ps-2] 

m 

[A ps“’^/^] 

K+ 

768.7 

152.5 

3.393 X 103 

7.094 X 10^ 

Na+ 

2.044 X 10^ 

166.1 

4.020 X 10^ 

1.717 X 10^ 

Ca2+ 

3.026 X 10^ 

190.2 

4.933 X 10^ 

1.874 X 10^ 

cr 

940.0 

189.7 

4.524 X 10^ 

1.061 X 10^ 


Table 3. Values of r)j, j = 1,2, 3, 4, estimated using all-atom MD simulations. 


estimated by calculating mean square displacements or velocity autocorrelation 
functions. In Table [21 we report the values of D which were estimated in [1 by 
calculating mean square displacements. 

Let us consider the coarse-grained model (|1.2I) - (|1.5I) and let (•) denotes an 
average over many realizations of a stochastic process. Multiplying equations (II.3p 
and (11.41) by Vi and C/j, respectively, we obtain the following ODEs for second 
moments: 

= 2{U,V), (3.3) 

= -2'ni{UiVi) + 2{UiZi). (3.4) 

Consequently, we obtain that {UiVi)=0 and {UiZi)=0 at steady state. 
Multiplying equations (ll.3p - (ll.5l) by V), Ui and Zj, and taking averages, we obtain 

= {uf)-r^i{Vl) + {VZ,), (3.5) 



Using {UiVi) = 0 and {UiZi) 


= {U,Zi) - 7?2(UZ,) - m{UiVi). (3.6) 


= 0, we obtain that {ViZi) 


{Uj 

{vd 


0 at steady state and 


(3.7) 


This equation is used in Table [3] to estimate rji using the MD averages (U?) and 
(V^) which are given in Table [2] Since we know the value of rji, we can also 
estimate the value of {Zf) by calculating the second moment of 


(Zf) 


fUi{t + At) 
V At 


Ui{t) 





(3.8) 


This value is reported in the last column of Table |21 Multiplying equation (II.4p 
by Zi and equation (|1.5p by Ui, we obtain 

^^{U^Z,) = (Zf) - m{V^Z,) - r]-,{u,z,) - Tjsiu!). 


(3.9) 
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Using {UiZi) = 0 and {ViZi) = 0, we obtain at steady state 


^ {ufy 

(3.10) 

Mnltiplying eqnation (11.21) hy Xi, Uj Ui and Zi and eqnations (ll.3l)-(ll.5l) by Xi 
and taking averages, we obtain the following system of ODEs for second moments: 

= 2(X,U,), 

(3.11) 

= {V^) + {X,U,), 

(3.12) 

at 

(3.13) 

^^{XiZi) = (UZi)-r?2(X,Zi)-r?3(X,C/i). 

(3.14) 

Conseqnently, we obtain at steady state (XjU) =D, {XiUi) = ■ 
rjiD and 

m{x^u,) V3{v^) 

{XiZ,) rjiD 


Using- (13.71) and (I3.10p. we have 


„ _{zl) ({v?)V 

(3.15) 

Finally, mnltiplying eqnation (11.51) by Zi, we obtain 


j^iZf) = -2rj2{Zf) - 2r?3(U,Z,) + r/i 

(3.16) 

Conseqnently, we obtain at steady state 


r]l = 2r]2{Zf). 


Using (I3.15p. we get 

y D {uf) • 

(3.17) 

lire valnes calcnlated bv (|3.10p. (|3.15p and (|3.17p are presented 

in Table [3l 


4. Accuracy of the coarse-grained model of ion 

The coarse-grained model (|1.2I) - (|1.5I) has fonr parameters rji, i = l,2,3,4. To 
parameterize this model, we have nsed fonr qnantities estimated from detailed 
MD simnlations, diffnsion constant D and steady state valnes of (V^^), {Uf) and 
(Zf). In particnlar, the coarse-grained model (ll.2p - (ll.5p will give the same valnes 
of these fonr qnantities, inclnding the valne of diffnsion constant D which is the 
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sole parameter of the BD model (11.11) . In this section, we explain why the coarse¬ 
grained description given by (m-dLS]) can be used as an intermediate model to 
couple BD and MD models. 


We begin by illustrating why Langevin dynamics (which is used in [l2[ for a 
similar multiscale problem) is not suitable for all-atom MD simulations studied 
in this paper. In |l2|, a few (heavy) particles with mass M and radius R are 
considered in the heat bath consisting of a large number of light point particles 
with masses m <C M. The collisions of particles are without friction, which 
means that post-collision velocities can be computed using the conservation of 
momentum and energy. In this case, it can be shown that the description of 
heavy particles converges in an apropriate limit to Brownian motion given by 
equation dni). One can also show that the model converges to Langevin dynamics 
(in the limit m/M —)• 0) [l^, [H! 1^ : 


dW = Vidt, 


(4.1) 


dVi = -'jVidt + dWi, i = 1,2,3, (4.2) 


where X = [Xi, X2, X3] is the position of a diffusing molecule, V = [Vi, V2, V3] is 
its velocity, D is the diffusion coefficient and 7 is the friction coefficient. In jl2|, 
Langevin dynamics dUll-dOl) is used as an intermediate model which enables 
the implementation of BD description (jl.ip and the original detailed model in 
different parts of the computational domain. 

Langevin dynamics (I0)-(I01) describes a diffusing particle in terms of its 
position and velocity, i.e. it uses the same independent variables for the description 
of an ion as the MD model dau-dMi). Langevin dynamics can be further 
reduced to BD model (|l.ip in the overdamped limit 7 —>^ 00 . However, it cannot 
be used as an intermediate model between BD and all-atom MD simulations 
considered in this paper, because it does not correctly describe the ion behaviour 
at times comparable to the MD timestep At. To illustrate this, let us parameterize 
Langevin dynamics (I4T])-(I01) using diffusion constant D and the second velocity 
moment {V^) estimated from all-atom MD simulations. To get the same second 
moment of velocity, Langevin dynamics requires that we choose 


7 = 



Discretizing equation (14.21) . the ion acceleration during one time step is 


(4.3) 


Vi{t + M)-Vi{t) 
At 


■'yVi{t) + 7 



(4.4) 


where [^ 1 , ^2 j Cs] is a vector of normally distributed random numbers with zero 
mean and unit variance. Using (|4.3p . the second moment of the right hand side 
of (14.41) is 


7 


2 



D2 D At 


(4.5) 


Using the MD values of D and (V)^) for K"*" which are given in Table |2] and 
using MD timestep At = 10“^ ps, we obtain that the second moment (|4.5p is 
equal to 4.44 x 10^ ps”"^. On the other hand, {Uf) estimated from all-atom 
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MD simulations and given in Table |2] is 4.86 x 10^ ps“^ which is one hundred 
times smaller. The main reason for this discrepancy is that Langevin dynamics 
postulates that the random force in equation (14.21) acting on the particle at time 
t is not correlated to the random force acting on the particle at time t + At. 
However, this is not true for all-atom MD simulations where random force terms 
at subsequent time steps are highly correlated. 

Since Langevin dynamics is not suitable for coupling MD and BD models, we 
need to introduce a stochastic model of ion behaviour which is more complicated 
than Langevin dynamics. The coarse-grained model (|1.2p - (|1.5p studied in this 
paper is a relatively simple example of such a model. Its parametrization, 
discussed in Section |3l guarantees that the coarse-grained model (iLl-dLS]) well 
approximates all-atom MD simulations at steady state. They both have the same 
value of diffusion constant D and steady state values of {V^), {Uf) and {Z‘f). Next, 
we show that the coarse-grained model (|1.2I) - (|1.5I) also compares well with all MD 
simulations at shorter timescales. We consider the rate of change of acceleration 
(jerk or the scaled derivative of force). We define the average jerk as a function 
of current velocity and acceleration of the ion: 


Jiv, u) = lim 

r —>-0 


{Ui{t + t) 


u \ Vi{t) = v,Ui{t) = u) 
r 


(4.6) 


To estimate J(u, u) from all-atom MD simulations, we calculate the rate of change 
of acceleration during each MD time step 


J{v, u) 


{Ui{t + At) - u I Vi{t) = V, Ui{t) = u) 
At 


(4.7) 


i.e. we run a long (nanosecond) MD simulation, calculate the values of {Ui{t + 
At) — Ui{t))/At during every time step and record their average in two-variable 
array J{v,u) indexed by binned values of Vi{t)=v and Ui{t) = u. Since the 
estimated J{v,u) only weakly depends on u, we visualize our results in Figured] 
using two functions of one variable, v, namely 


= J(^;,0), 


and 


Mv) 


'CO 

J{v,u)pu{u) du, 

J —oo 


(4.8) 


where Pu{u) is the steady state distribution of Ui estimated from the same 
long time MD trajectory. As before, we use all three dimensions to calculate 
the averages J{v,u) and Pu{u)- Function Ji{v) (which gives jerk at the most 
likely value of Ui) is plotted using crosses and function J 2 {v), the average over 
Ui variable, is plotted using circles in Figure dl In order to compare all-atom 
MD simulations with the coarse-grained model (fL^-lOl). we calculate the 
corresponding jerk matrix J{v,u) for the coarse-grained model. We denote by 
p{v,u,z) the stationary distribution of the stochastic process (|1.3l) - (|1.5p . i.e. 
p{v, u, z) du dudz is the probability that Vi{t) G [u, u -|- du), Ui{t) G [u, u -|- du) and 
Zi{t) G [z, z -|- du). Then the jerk matrix (14.61) of the coarse-grained model (11.21) - 
(frb is 


J(u, u) = 


lim 

r —>-0 


{Ui{t + t) -u \ Vi{t) = V, Ui{t) = u, Zi{t) = z) 


p{v, u, z) dz. 


— OO 


r 
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Figure 1. Comparison of the rate of change of acceleration (jerk) computed by all-atom MD 
simulations and by the coarse-grained model (fT:B-(IT31). MD results are visualized using 
functions Ji(i>) (black crosses) and J 2 {v) (blue circles) defined by equation (14.81) . The result 
obtained by the coarse-grained model is given by formula (14.101) (red solid line). We consider 
(a) ion; (b) Na"*" ion; (c) Ca^'*" ion and (d) Cl“ ion. Parameters are given in Tahles^and^ 


Using (11.41) , we rewrite it as 

'OO 

J{v,u)= {—r]iv + z)p{v,u,z)dz. 


(4.9) 


The stationary distribntion p{v, u, z) of (jl.3p - (|1.5l) is Ganssian with mean [0,0, 0]^ 
and stationary covariance matrix: 


1 


vl 0 


0 


0 pi'ni 0 


2??ir?2r?3 \ 0 0 

Conseqnently, eqnation (|4.9p implies 

J(n, u) = —rjiv. 


vimvL 


(4.10) 
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ion 

Ai 

[ps-i] 

A2 

[ps-i] 

A3 

[ps-i] 

q 

[ps] 

t*2 

[ps] 

K+ 

-127.0 

-12.75-b 27.58 i 

-12.75 - 27.58 i 

3.08x10-2 

-9.39x10-3 

Na+ 

-140.1 

-12.99-b 47.47 i 

-12.99-47.47i 

6.15x10-3 

-2.35x10-2 

Ca 2 + 

-163.1 

-13.58-b 57.84 i 

-13.58 - 57.84 i 

1.47x10-3 

-2.48x10-2 

cr 

-162.9 

-13.41 -b 30.25 i 

-13.41 - 30.25 i 

2.50x10-2 

-1.13x10-2 


Table 4. Eigenvalues Xj, j = 1, 2, 3, of matrix B defined by (1^ and time shifts t* and tl- 
Symbol i denotes the imaginary unit. 


In Figure [U we plot (14.101) using the red solid line. The comparison with all atom 
MD results (circles and squares) is excellent for all four ions considered in this 
paper. In particular, we have shown that the coarse-grained model dLl-dlSl) 
provides a good description of the rate of change of acceleration (jerk) at the MD 
timescale. We make use of this property in Section [ 6 ] where we use the same time 
step (At = 10“^ ps) for both the coarse-grained model (|1.2p - (|1.5l) and all-atom 
MD simulations. The coarse-grained model (ll.2l) - (|1.5p can also be coupled with 
BD description (|l.ip . which uses much larger time steps, as we show in the next 
section. 


5. Prom the coarse-grained model ()1.2p - (|1.5p to Brownian dynamics 


Let us consider the three-variable subsystem dLSlI-dlSl) of the coarse-grained 
model. Denoting y* = [V),17i,Zj], equations (|1.3p - (|1.5l) can be written in vector 
notation as follows 

dyi = ByiAt + hAWi, (5.1) 

where matrix B E and vector b E are given as 



(5.2) 


Let us denote the eigenvalues and eigenvectors of B as Xj and Vj = [I'lj, 1 ^ 2 j, 1 ^ 3 j]: 
j = 1,2,3, respectively. The eigenvalues of B are the solutions of the characteristic 
polynomial 

+ r]2X^ + {rji + r]s)X + mm = 0 . 


Since 771 , m and r /3 are positive parameters, we conclude that real parts of all 
three eigenvalues are negative and lie in interval (—172,0). Using the values of r]j, 
j = 1,2,3, given in Table [3l we present the values of eigenvalues Xj, j = 1,2,3, in 
Table m The eigenvalues Xj, j = 1,2,3, are distinct. The general solution of the 
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SDE system (15.11) can be written as follows 


yi{t) = ^{t)c + ^{t) 


't 

0 


<^-\s)hdWi, 


(5.3) 


where c E is a constant vector determined by initial conditions and matrix 
is given as <l>(t) = [exp(Ait)^'i | exp(A 2 t )^'2 | exp(A 3 t)i/ 3 ], i.e. each 
column is a solution of the ODE system dyj = Byi dt. Considering deterministic 
initial conditions, equation (15.31) implies that the process is Gaussian at any time 
t > 0. Equations for means, variances and covariances then uniquely determine the 
distribution of yi{t) for t > 0. Equations for means can be written in the vector 
form as d(yj) = B (yj) dt. Equations for variances and covariances are given in 
Section |3] as equations (I3.3p - (l3.6p . (13.9p . (I3.1ip - (l3.14p and (I3.16p . 

There are two important conclusions of the above analysis. First of all, 
eigenvalues Xj, j = 1,2,3, given in Table |4] satisfy 


Ai < Re A 2 — ReA 3 < 0, 


where Re denotes the real part of a complex number. There is a spectral gap 
between the hrst eigenvalue and the complex conjugate pair of eigenvalues. If we 
used this spectral gap, we could reduce the system to two evolution equations 
for times t;:|>l/|Ai|. However, there is no spectral gap to reduce the system to 
Langevin dynamics (|4.1I) - (|4.2I) . In particular, we again confirm our conclusion 
that a coarse-grained approximation of ion behaviour is not given in terms of 
Langevin dynamics. Our second conclusion is that on a picosecond time scale, we 
can assume stationarity in (15.11) to get 

dXi = -^dWi, i = 1,2,3. (5.4) 

mm 

Using (13.7p . (|3.15l) and (13.171) . we have 

mm 

Consequently, equation (15.41) is equivalent to BD description (II.ip . The 
convergence of (|1.2I) - (|1.5I) to the BD model is illustrated in Figure [2j a). We solve 
the system of 10 ODEs for variances and covariances given as equations dSS])- 
(j3.6l) . (|3.9p . (|3.1ip - ()3.14l) and (|3.16l) . We consider (deterministic) zero initial 
conditions, i.e. Xj(0) = Vi{d) = 7/^(0) = Zi{ll) = 0. All moments are then initially 
equal to zero. We plot the mean square displacement {Xf) as a function of time. 
We compare it with the mean square displacement of BD model (|l.ll) which is 
given as 2Dt. We observe that there is an approximately constant shift, denoted t^, 
between both solutions for times t > 0.2 ps. We illustrate this further by plotting 
{Xf{t + tj)) in Figure [2Ka). The values of shift for different ions estimated by 
solving the ODEs for second moments with zero initial conditions are given in 
Table HI 

Next, we show how the BD model dni) and the coarse-grained model m- 
(jl.5l) can be used in different parts of the computational domain. This coupling 
will form one component of multiscale methodology developed in Section |6l 
BD algorithms based on equation dni have been implemented in a number 
of methods designed for spatio-temporal modelling of intracellular processes. 
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time t [ps] 


(b) 



Xi [A] 


Figure 2. (a) Comparison of the coarse-grained model (n31)-(fT31l and BD description dH) for 
K+ ion. The mean square displacement computed by solving 10 DDEs (13.31) - (13.6L (13. 9L (13.111) - 
(isna and (I3A61) with zero initial conditions (black solid line). The same curve shifted by the 
value of t) is plotted as a red dashed line, (b) Test of accuracy of the multiscale approach 
in fla U fl 4 U Sis for K'*' ion. Histogram obtained by simulating 10® ions over time 10® ps is 
compared with the analytical result (15.61) (red solid line). Grey bars show the ion density in SI 3 , 
the green bar shows the ion density in SI 4 and blue bars show the ion density in SI 5 . Parameters 
are given in Tahles^and^ 


including Smoldyn j28|, MCell j29|| and Green’s-function reaction dynamics 
Smoldyn discretizes dni) using a fixed BD time step AT, i.e. it computes the 
time evolution of the position X = X(t) of each molecule by 


Xi{t + At) = Xi{t) + V2DAT(i, i = 1,2,3, 


(5.5) 


where [^ 1 , ^2 j Cs] is a vector of normally distributed random numbers with zero 
mean and unit variance. We use discretization (|5.5I) of BD model (|l.ip in this 
paper. BD time step AT has to be chosen much larger than the MD time step 
At. We use AT = 0.5 ps, but any larger time step would also work well. We could 
also use a variable time step, as implemented in the Green’s Function Reaction 
Dynamics [30l |. 

In Section | 6 l we consider all-atom MD simulations in domain Our 

main goal is to design a multiscale approach which can compute spatio-temporal 
statistics with the MD-level of detail in relatively small subdomain Di C D by 
using BD model (15.5p in the most of the rest of the computational domain. This 
is achieved by decomposing domain D into hve subdomains Qj, j = 1,2,3,4, 5 
(see equation m and discussion in Section El). We use MD in Di, the coarse¬ 
grained model dH-dls]) in D 3 and the BD model (|5.5I) in D 5 . The remaining 
two subdomains, D 2 and D 4 , are two overlap (hand-shaking) reg ions where two 
different simulation approaches can be used at the same time |12 . kMI | . In the rest 
of this section, we focus on simulations in region D3 U D4 U D5 which concerns 
coupling the coarse-grained model (ll.2p - (ll.5l) with the BD model (15.51) . We use the 
coarse-grained model in D 3 U D4 and the BD model (15.51) in D4 U D 5 . In particular, 
we use both models in the overlap region D 4 . Each particle which is initially in 
D 3 is simulated according to (|1.2p - (|1.5p (discretized using time step At) until it 
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enters ^ 5 . Then we nse (15.51) to evolve the position of a particle (over BD time 
steps of length AT) nntil it again enters ^3 when we switch the description back 
from the BD model to the coarse-grained model. In order to do this, we have to 
initialize variables V), Ui and Zj, i = 1, 2,3. We nse deterministic initial conditions, 
Vi(0) = C/i(0) = Zj(0) = 0, disccnsed above. 

In Fignre Efb), we present an illnstrative simnlation where D 3 U D 4 U D 5 = 
for simplicity. We nse D 3 = (h, 00) x M^, D 4 = [—/i, h\ x and D 5 = (—00, —h) x 
M^, where h = lK. We report averages over 10® simnlations of ions, half of 
them are initiated at X(0) = [/i, 0,0], i.e. they initially follow the coarse-grained 
model (IL21)-(II3|) with zero initial condition for other variables (Vi(0) = Ui{Q) = 
Zj(0)=0). The second half of ions are initiated at X(0) = [—/i, 0, 0], i.e. they 
initially follow BD description (|5.5I) . We plot the (marginal) distribntion of ions 
along the first coordinate (Xi) at time 10^ ps in Fignre E^b). The compnted 
histogram is plotted nsing bins of length 2 A, i.e. the overlap region D 4 is eqnal 
to one bin (visnalized as a green bar). Grey (resp. bine) bars show the density 
of ions in D 3 (resp. D 5 ). We compare onr resnlts with the analytical distribntion 
compnted for BD description (11.11) at time t = 10^ ps given by 


q{xi) = 


10 ® 


A^/irDt 


f 

{xi-hf 


(xi -k h)^ 

(exp 

4Dt 

+ exp 

4Dt 


(5.6) 


The compnted histogram compares well with (15.6p . althongh we can observe a 
small error: the green bar is slightly taller than the corresponding valne of (15.61) . 
If we wanted to fnrther improve the accnracy, we conld take into acconnt that 
there is time shift discnssed above, introdnced to the mnltiscale approach 
by nsing the deterministic initial conditions, V)( 0 ) = C/i( 0 ) = Zj( 0 ) = 0 , for ions 
entering domain D 3 . Another possibility is to sample the initial condition for V), 
Ui and Zi from a snitable distribntion. If we nse the stationary distribntion of 
snbsystem (fT3])-(0]). then iy^) does not evolve and is eqnal to 


{V?) = 


vl 




Snbstitnting this constant for (V^) into (13.121) . the system of 10 ODEs for second 
moments of (|1.2p - (|1.5l) simplifies to 4 ODEs (I3.1ip - (l3.14p . Solving system (13.111) - 
(|3.14l) with zero initial conditions (assnming Xi (0) = 0), we can again compnte 
the mean sqnare displacement. As in Fignre EJa), it can be shifted in time to 
better match with the BD resnlt, 2Dt. We denote this time shift as Its valnes 
are given in Table |4] We observe that is negative and is positive for all fonr 
ions considered in Table E] Both time shifts and (together with optimizing 
size h of the overlap region) conld be nsed to fnrther improve the accnracy of 
mnltiscale simnlations in D 3 U D 4 U D 5 jl2|. However, onr main goal is to introdnce 
a mnltiscale approach which can nse all-atom MD simnlations in Di. Since MD 
simnlations are compntationally intensive, we will only consider 100 realizations 
of the mnltiscale method in Section [ 6 l In particnlar, the Monte Carlo error will 
be larger than the error observed in Fignre Elb)- Thns, we can nse the above 
approach in D 3 U D 4 U D 5 withont introdncing observable errors in the mnltiscale 
method developed in the next section. 
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^5 


^3 



Ion simulated using all-atom MD 
in Hi. 

Ion described by the coarse-grained 
model (ll.2l) - (ll.5l) in 112 ) ^^3 and Il 4 . 

Ion described by BD model (|5.5p 
in Il 4 and Us. 

If an ion is in Hi or H 2 , then water 
molecules are simulated using 
all-atom MD in Hi. 


Figure 3. Schematic of multiscale set up. Note that the schematic is drawn in two spatial 
dimensions to enable better visualization, but all models are formulated and simulated in three 
spatial dimensions. 


6. Coupling all-atom MD and BD 

Let us consider all-atom MD in domain H C which is so large that direct 
MD simulations would be too computationally expensive. Let us assume that 
a modeller only needs to consider the MD-level of detail in a relatively small 
subdomain Hi C H while, in the rest of the computational domain, ions are 
transported by diffusion and BD description is applicable. For example, 

domain Hi could include binding sites for ions or (parts of) ion channels. In 
this paper, we do not focus on a specific application. Our goal is to show that 
the coarse-grained model dOll-dLSI) is an intermediate model between all-atom 
MD and BD which enables the use of both methods during the same dynamic 
simulation. To achieve this, we decompose domain H into five subdomains, 
denoted Hj, j = l,2,3,4,5, see Figure El These sets are considered pairwise 
disjoint (i.e. Hj n Hj = 0 for j) and 

H = Hi U H2 U H3 U H4 U H5. ( 6 . 1 ) 

In our illustrative simulations, we consider the behaviour of one ion. If the ion is 
in subdomain Hi, then we use all-atom MD simulations as described in Section [2] 
In particular, the force between the ion and a water molecule is obtained by 
differentiating potential (12.11) . provided that the distance between the ion and 
the water molecule is less than the cutoff distance (T/2). Let us denote the force 
exerted by the ion on the water molecule by Fj^(rjo, r'ji, rj 2 ), where r^o (resp., 
rn and ri 2 ) is the distance between the ion and the oxygen site (resp., the first 
and second hydrogen sites) on the water molecule. We use periodic boundary 
conditions for water molecules in Hi. 

Whenever the ion leaves Hi, it enters H 2 where we simulate its behaviour using 
the coarse-grained model (UlIl-dLSl). We still simulate water molecules in Hi and 
we allow them to experience additional forces exerted by the ion which is present 
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in ^ 2 - These forces have the same functional form, Fj^, as in MD, but they have 
modified arguments as follows 

Fj«,(rjo + wdist(X, ni),rji + w dist(X, fli), rj 2 + a;dist(X, fli)), (6.2) 

where a; > 0 is a parameter and dist(X, fli) is the (closest) distance between 
the ion at position X and subdomain fli. If the ion is in region fls U U fls, 
then water molecules in fli are no longer simulated. We use the coarse-grained 
model (|1.2p - (|1.5p to simulate the ion behaviour in and the BD model (15.51) in 
fis. Overlap region 124 is used to couple these simulation methods as explained in 
Section |5l 

In Section |5l we have already presented illustrative simulations to validate the 
multiscale modelling strategy chosen in region Us U 04 U Us. Next, we focus on 
testing and explaining the multiscale approach chosen to couple region 121 with 
122 - The key idea is given by force term (16.21) which is used for MD simulations of 
water molecules in I2i when an ion is in 122 . This force term has two important 
properties: 


(i) If an ion is on the boundary of I2i, i.e. XG(9I2 i, then dist(X, I2i) =0 
and force (| 6 . 2 I) is equal to force term Fj^(rjo,r*!,^* 2 ) used in I 2 i. 

(ii) If a;dist(X, I2i) > L/2, then force (| 6 . 2 p is equal to zero. 

Property (i) implies that formula (16.21) continuously extends the force term used 
in MD. In particular, water molecules do not experience abrupt changes of forces 
when the ion crosses boundary (9I2i. Property (ii) is a consequence of the cutoff 
distance used (together with the reaction field correction Q) to treat long-range 
interactions. In our illustrative simulations, we use 


Di — 


L L 

3 

■ L L L L ■ 

~2' 2 

and 122 = 

~2 ~ ^'2 ^ 


\I2i. 


(6.3) 


Property (ii) implies that extra force (16.21) is equal to zero on boundary 5122 \ clI2i 
which is the boundary between regions 122 and This is consistent with the 
assumption that ions in region Ds U 124 U ^^5 do not interact with water molecules 
in region 121 . 

If an ion is in I2i, we use all-atom MD as formulated in Section |2] Periodic 
boundary conditions are implemented in MD simulations. Water molecules are 
subject to forces exerted by the ion at its real position X in I2i, but also by 
its copies at periodic locations X-|- {iL,jL,kL) where i,j,k^Z. When the ion 
moves to 122 , one of its copies is in I2i. Force term (16.2p is designed in such a way, 
that the strength of interaction decreases (for every copy of the ion) with the 
distance, dist(X, I2i), between the real position of the ion and I2i. In particular, 
force term (j 6 . 2 l) ensures that there are continuous changes of all forces when the 
ion moves between regions 121 , 122 and ^ 3 . 

In Figure SJ we present results of simulations of K"*" ion in region I2i U 122 . 
We consider 100 realizations of a multiscale simulation with one ion. Its initial 
position is X(0) = [—L/2, 0,0] which lies on boundary clI2i. We simulate each 
realization for time 10 ps which is short enough that all trajectories stay inside 
the ball of radius L/2 centred at X(0). Then Xi-coordinate of the trajectory 
determines whether the ion is in I 2 i or 122 . If f^i(f) ^ ~T/ 2 , then the ion is in 
I2i and it is simulated using all-atom MD. If X 2 {t) < —L/2, then the ion is in 122 
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(a) 



0 2 4 6 8 10 

time [ps] 


(b) 



Figure 4. (a) One hundred realizations of a multiscale simulation of ion initiated at 
[—L/2, 0,0]. We plot Xi coordinate as a function of time. Ion is described by all-atom MD for 
X\ > —1//2 and by the coarse-grained model (II.2ll - (11.51) for Xi < —L/2. The boundary between 
Hi and 0.2 is visualized using the black dashed line. We use tj = 1 ire (16.21) . (b) The mean sguare 
displacement in the first coordinate of ion simulated in Hi U O 2 and computed as the average 
of 100 realizations for tj = 1 (blue circles), uj = 2 (black crosses) and tj = 10 (green sguares). 


and evolves according to the coarse-grained model (ll.2l) - (ll.5l) . In Figure Ufa), we 
use (|6.2p with w = 1 and plot Xi coordinates of all 100 realizations. We observe 
that the computed trajectories spread on both sides of boundary Srii (dashed 
line) without any significant bias. The mean square displacement is presented in 
Figure SKb) for three different values of uj. The results compare well with 
which is the mean square displacement of one coordinate of the diffusion process. 

We conclude with illustrative simulations which are coupling all-atom MD with 
BD. We use domain D G decomposed into five regions as in equation (16.11) . 
where Di and D 2 are given by (16.3p . and 


D3 

D4 


L L 


L L 


-— — — h\ , — -)- — -)- hi 

2 2a; ’ 2 2a; 


\ (Di U D 2 ) 
1 3 


L L , L L , 

-— — — — '^ 2 1 — “t“ — H” h} /^'2 

2 2c<j 2 2c<; 


\ (Di U D 2 U D 3 ), 


(6.4) 

(6.5) 


D 5 — \ (D 3 U D 2 U D 3 U D 4 ), (b' 6 ) 

where a; = 10, /ii = L/20 and /i2=T/10. Then the BD domain is D 5 =]R^\ 
[—7L/10, 7L/10]^. We place an ion at the origin (centre of MD domain Di), 
i.e. X(0) = [0,0,0], and we simulate each trajectory until it reaches the distance 
4L = 99.32 A from the origin. Let T{r) be the time when a trajectory first reaches 
distance r from the origin. In Figure |5l we plot escape time T{r) as a function 
of distance r. We plot the value of T{r) for each realization as a blue point. The 
largest computed escape times (for r = AL) are 38, 506ps for K"*" and 47, 212ps 
for Na"*". They are outside the range of panels in Figure |5l but the majority of 
data poins are included in this figure. We also plot average {T{r)) (red solid 
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(a) (b) 




distance [A] distance [A] 

Figure 5. Escape time T{r) to reaeh distance r from the origin computed by the multiscale 
method. We consider (a) K’^ ion; and (b) Na"*" ion. We plot escape times for individual 
realizations (blue points), the mean escape time estimated from 100 realizations (red solid line) 
and the theoretieal 95% confidence interval (16.711 (green area). We use a; = 10 in (16.211 . 


line) together with 95% confidence intervals. They are compared with theoretical 
results obtained for the BD model m- The escape time distribution for the 
BD model dni has mean equal to i(T{r)) = L?/{QD) and standard deviation 
L^/(3 \/ToZ 1). The corresponding theoretical 95% confidence interval (for 100 
samples) is 


L2 

W - 


L2 \ 

^ ^’^^300 J ■ 


(6.7) 


This interval is visualized as the green area in Figure [S] We note that it would 
be relatively straightforward to continue the presented multiscale computation 
and simulate ion diffusion in domains covering the whole cell. The most 
computationally intensive part is all-atom MD simulation in Di U D 2 . However, 
once the ion enters H 5 , we can compute its trajectory very efficiently. We could 
further increase the BD time step in parts of D 5 which are far away from or 
we could use event-based algorithms, like Green’s-function reaction dynamics j30| 
or First-passage kinetic Monte Carlo method |32l | . to compute the ion trajectory 
in region D 5 . 


7. Discussion 

In this paper, we have introduced and studied the coarse-grained model m- 
([LSl) of an ion in aqueous solution. We have parameterized this model using 
all-atom MD simulations for four ions (K+, Na"*", Ca^'*' and Cl“) and showed that 
this model provides an intermediate description between all-atom MD and BD 
simulations. It can be used both with MD time step At (to couple it with all-atom 
MD simulations) and BD time step AT (to couple it with BD description (ll.ip i. 
In particular, the coarse-grained model enables multiscale simulations which use 
all-atom MD and BD in different parts of the computational domain. 
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In Section | 6 l we have illustrated this multiscale methodology using a first 
passage type problem where we have reported the time taken by an ion to reach 
a specific distance. Possible applications of this multiscale methodogy include 
problems where a modeller considers all-atom MD in several different parts of 
the cell (for example, close to binding sites or ion channels) and wants to use 
efficient BD simulations to transport ions by diffusion between regions where MD 
is used. The proposed approach thus enables the inclusion of MD-level of detail 
in computational domains which are much larger than would be possible to study 
by direct MD simulations. 

Although the illustrative simulations in Section | 6 ] are reported over distances 
of the order of 10 ^ A, this is not a restriction of the method. Most of the 
computational time is spent by considering all-atom MD in Di U D 2 . BD uses 
much larger time step which enables us to father extend BD region D 5 (and 
consequently, the original domain D). Moreover, if we are far away from MD 
domain Di, we can further increase the efficiency of BD simulations b y u sing 


different BD time steps in different parts of the BD subdomain D 5 jl2|, or 
by using event-based BD algorithms j^, |32|. The computational intensity of 
BD simulations can be further decreased by using multiscale methods which 
efficiently and accurately combine BD models with lattice-based (compartment- 
based) models (^ . . Such a strategy have been previously used for modelling 

intracellular calcium dynamics P, U or actin dynamics in filopodia (35| . and 
enables us to extend both temporal and spatial extent of the simulation. 

In the literature, MD methods have been used to estimate parameters of BD 
simulations of ions (3^ . There has also been a lot of progress in systematic coarse- 
graining of MD simulations The approach presented in this paper not only 
uses all-atom MD simulations to estimate parameters of a coarser description, 
but it also designs a multiscale approach where both methods are used during 
the same simulation. Methods which adaptively change the resolution of MD on 
demand have been previously reported in (^, . They include algorithms which 

couple all-atom MD with coarse-grained MD. The coarse-grained model developed 
in this work does not include any water molecules and has different application 
areas. One of them is modelling of calcium induced calcium release through IP 3 R 
channels [3| which is discussed as a motivating example in Introduction. MD 
simulations in this paper use the three-site SPC/E model of water. An open 
question is to extend our observations and analysis to other MD models of water, 
which include both more detailed water models with additional sites |lfil . Il7l | and 
coarse-grained MD models of water 
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